
afripolis_shp = file.path(wkdir,'africapolis')

require(Rcpp)
setwd(wkdir)
output_version <- '2020-09-02'
 

# Maiduguri coordinates
africapolis_dir = paste0(wkdir,'/local/gis_data/003_boundaries')
maiduguri <- st_read(dsn=paste0(africapolis_dir,'/boko_haram_capital'),layer='Maiduguri_africapolis_pt')
bhstate <- st_read(dsn=paste0(wkdir,'/local/gis_data/003_boundaries'),layer='bhstate_ln')

## africapolis, poly
## https://www.africapolis.org/data
afripolis <- st_read(dsn=dirname(afripolis_shp),layer=basename(afripolis_shp))
afripolis <- subset(afripolis,ISO %in% c('CMR','NER','NGA','TCD'))
afripolis_pt = st_centroid(subset(afripolis,select=c("pop2015","pop1950","pop1960","pop1970","pop1980","pop1990","pop2000","pop2010")))
names(afripolis_pt)

# select capitals from afripolis data #
capitals = c('Abuja','Niamey','Ndjamena/Kouss?ri [TCD]','Yaound?')
cap = subset(afripolis,Name %in% capitals)
cap_pt = st_centroid(cap)

# select primary cities from afripolis data #
pri = c('Lagos','Niamey','Ndjamena/Kouss?ri [TCD]','Yaound?')
pri = subset(afripolis,Name %in% pri)
pri_pt = st_centroid(pri)

port = c('Douala')
port = subset(afripolis,Name %in% port)
port_pt = st_centroid(port)

# https://www.prio.org/data/11
petro_shp = file.path(petro_dir,"PETRO_Onshore_080907")
petro_ply = st_read(dsn=dirname(petro_shp), layer=basename(petro_shp),stringsAsFactors = FALSE)
oil_ply = subset(petro_ply,COUNTRY %in% c('Cameroon','Chad','Niger','Nigeria'))

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
if ((length(unique(ntl.zones$OBJECTID))==dim(ntl.zones)[1])==FALSE) { break }

##
## join pop data
##
require(dplyr)
pop.zones <- st_join(ntl.zones, afripolis_pt, join = st_contains)  %>% 
  group_by(OBJECTID) %>%
  summarize(pop1950 = sum(pop1950),pop1960=sum(pop1960),pop1970=sum(pop1970),pop1980=sum(pop1980),pop1990=sum(pop1990),pop2000=sum(pop2000),pop2010=sum(pop2010),pop2015=sum(pop2015))

# add area km2 #
ntl.zones$sf_area_km2 <- st_area(ntl.zones) / 1000000


## poly to center point ##
albers_africa = "+proj=aea +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"
ntl.zones.prj <- ntl.zones %>%
  # 102022 = WGS 1984 Albers for Africa
  st_transform(albers_africa) %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)

## 3 definitions of boko haram area/polygon
bh3s = read_sf(dsn=paste0(wkdir,"/local/gis_data/003_boundaries"),layer="bh3s")
bhbr = read_sf(dsn=paste0(wkdir,"/local/gis_data/003_boundaries"),layer="bhbr_dst")
bhb = read_sf(dsn=paste0(wkdir,"/local/gis_data/003_boundaries"),layer="bhb")

ntl3.dt <- dplyr::bind_cols(ntl.zones,ntl.zones.pts)

(require(nngeo))

  # dist2maiduguri : boko haram capital
  ntl3.dt$dist2bhc <- unlist(nngeo::st_nn(ntl.zones.pts, maiduguri, k = 1, returnDist = T)$dist) / 1000
  summary(ntl3.dt$dist2bhc)
  
  ntl3.dt$negdist2bhc = -1 * ntl3.dt$dist2bhc
  
  #dist2bordbh3s boko-haram 3 states
  ntl3.dt$dist2bordbh3s <- unlist(nngeo::st_nn(ntl.zones.pts, bh3s, k = 1, returnDist = T)$dist) / 1000
  
  #dist2bordbhbr boko-haram between rivers
  ntl3.dt$dist2bordbhbr <- unlist(nngeo::st_nn(ntl.zones.pts, bhbr, k = 1, returnDist = T)$dist) / 1000
  
  #dist2bordbhb boko-haram borno only
  ntl3.dt$dist2bordbhb <- unlist(nngeo::st_nn(ntl.zones.pts, bhb, k = 1, returnDist = T)$dist) / 1000
  
  # dist2cap
  ntl3.dt$dist2cap_cmr <- unlist(nngeo::st_nn(ntl.zones.pts, cap_pt[cap_pt$ISO=='CMR',], k = 1, returnDist = T)$dist) / 1000
  
  # dist2pri
  ntl3.dt$dist2pri_cmr <- unlist(nngeo::st_nn(ntl.zones.pts, pri_pt[pri_pt$ISO=='CMR',], k = 1, returnDist = T)$dist) / 1000

  # dist2port
  ntl3.dt$dist2port_cmr <- unlist(nngeo::st_nn(ntl.zones.pts, port_pt[port_pt$ISO=='CMR',], k = 1, returnDist = T)$dist) / 1000

dist.dt <- as.data.frame(ntl3.dt)

dist.dt2 <- dist.dt %>%
  dplyr::select("OBJECTID","GID_0","dist2cap_cmr","dist2pri_cmr","dist2port_cmr")

save.dta13(dist.dt2, 
           paste0(wkdir,"/proc_data/DIST_cmr",output_version,".dta"))